O’Hara on GitHub


knitr::opts_chunk$set(echo = TRUE, message = FALSE, warning = FALSE)

library(raster)
library(oharac)
oharac::setup()
source(here('fb_slb_fxns.R'))

1 Summary

Combine vulnerability scores with AquaMaps species information and functional groups.

2 Data

See individual trait scripts.

3 Methods

3.1 Gather and assemble data

spp_fxn_groups <- read_csv(here('_data/functional_groups_from_traits.csv'))

spp_info <- get_am_spp_info()

spp_by_gp <- spp_info %>% 
  dplyr::select(am_sid, sciname) %>%
  inner_join(spp_fxn_groups, by = c('sciname' = 'species')) ### 18285 spp matched in fxn gps

spp_vuln <- get_spp_vuln() %>%
  select(sciname = species, stressor, score, sd_score) %>%
  distinct()

spp_vuln_am <- spp_info %>%
  inner_join(spp_vuln, by = c('sciname'))
nspp <- spp_vuln_am$sciname %>% n_distinct() 
### 10016 species with vulnerability scores in the AquaMaps dataset overall
### (will improve when AquaMaps names are rectified with WoRMS)

spp_vuln_by_gp <- spp_by_gp %>%
  inner_join(spp_vuln, by = c('sciname'))

nspp <- spp_vuln_by_gp$sciname %>% n_distinct() 
### 7985 species with both functional groups and vulnerability scores.  Of
### these, the highest group-gapfill level is ~3.2, not so bad...

3.2 Examine vulnerability by stressor and functional group

if(file.exists(here('figs/vuln_plot_by_group.png'))) {
  message('plot exists, not going to regenerate it!')
} else {
  spp_vuln_variation <- spp_vuln_by_gp %>%
    group_by(stressor) %>%
    mutate(diff = score - mean(score)) %>%
    group_by(p_group) %>%
    mutate(gp_lbl = paste0(p_group, ' (n = ', n_distinct(sciname), ')')) %>%
    ungroup()
  spp_vuln_var_summary <- spp_vuln_variation %>%
    group_by(stressor, gp_lbl) %>%
    summarize(mean = mean(diff),
              sd   = sd(diff)) %>%
    ungroup()
  x <- ggplot() +
    theme_minimal() +
    theme(axis.text = element_text(size = 7)) +
    geom_jitter(data = spp_vuln_variation, 
                aes(x = stressor, y = diff), color = 'grey70', alpha = .1) +
    geom_hline(yintercept = 0, color = 'red', size = 1) +
    geom_linerange(data = spp_vuln_var_summary, 
                   aes(x = stressor, ymin = mean - sd, ymax = mean + sd), 
                   color = 'black', size = .5) +
    geom_point(data = spp_vuln_var_summary, shape = 21, size = 1.5,
               aes(x = stressor, y = mean), color = 'black', fill = 'yellow') +
    geom_hline(yintercept = 0, color = 'red', size = .25) +
    coord_flip() +
    theme(axis.title.y = element_blank()) +
    facet_wrap(~gp_lbl) +
    labs(y = 'difference (vuln score - mean(vuln score))')
  
  ggsave(here('figs/vuln_plot_by_group.png'), height = 8, width = 8)
}

3.3 Map CHI for one species

Select a species (or a small handful) and calculate cumulative impacts across range.

  • Dermochelys coriacea (leatherback turtle): Rep-3437
  • Poroderma pantherinum (Leopard catshark): Fis-23221
  • Megaptera novaeangliae (humpback whale): ITS-Mam-180530
  • Thunnus orientalis (Pacific bluefin tuna): Fis-123736
  • Clupea harengus (Atlantic herring): Fis-29344

While we have stressor data from BD/CHI at 10km x 10km resolution, species distributions are at 0.5° resolution… for now let’s use some test rasters (from setup_stressors.Rmd) of the stressors in their most recent available year.

3.3.1 The stressors

str_maps <- data.frame(f = list.files(here('_spatial/stressors_hcaf'), full.names = TRUE)) %>%
  mutate(str_chi = str_remove(basename(f), '_[0-9]{4}.+'),
         year = str_extract(basename(f), '[0-9]{4}') %>% as.integer()) %>% 
  full_join(read_csv(here('_raw/strs_vuln_to_chi.csv')) %>%
              mutate(str_chi = str_split(str_chi, ';')) %>%
              unnest(str_chi)) %>%
  filter(!is.na(str_vuln) & !is.na(str_chi) & str_vuln != 'biomass_removal')

str_stack <- raster::stack(str_maps$f) %>%
  setNames(str_maps$str_vuln)
### fill NAs with zero for math reasons
values(str_stack)[is.na(values(str_stack))] <- 0

plot(str_stack)

spp_vec <- c('Dermochelys coriacea (leatherback turtle)'  = 'Rep-3437', 
             'Carcharodon carcharias (great white shark)' = 'Fis-23071', 
             'Megaptera novaeangliae (humpback whale)'    = 'ITS-Mam-180530', 
             'Thunnus orientalis (Pacific bluefin tuna)'  = 'Fis-123736', 
             'Clupea harengus (Atlantic herring)'         = 'Fis-29344')

test_vuln_df <- spp_vuln_by_gp %>%
  filter(am_sid %in% spp_vec) %>%
  select(am_sid, sciname, p_group, stressor, score, sd_score) %>%
  filter(stressor %in% str_maps$str_vuln)

test_maps_df <- get_am_spp_cells(prob = 0) %>%
  filter(am_sid %in% spp_vec) %>%
  mutate(presence = 1)

for(i in seq_along(spp_vec)) {
  ### i <- 2
  spp <- spp_vec[i]
  
  message('processing impacts for ', names(spp_vec)[i])
  cat(sprintf('\n\n### %s \n', names(spp_vec)[i]))
  tmp_spp_range_rast <- test_maps_df %>%
    filter(am_sid == spp) %>%
    map_to_hcaf(which = 'presence')
  tmp_spp_vuln_df <- test_vuln_df %>%
    filter(am_sid == spp) %>%
    filter(score > 0)
  
  impact_list <- vector('list', length = nrow(tmp_spp_vuln_df)) %>%
    setNames(tmp_spp_vuln_df$stressor)
  for(j in 1:nrow(tmp_spp_vuln_df)) {
    ### j <- 1
    s <- tmp_spp_vuln_df$stressor[j]
    v <- tmp_spp_vuln_df$score[j]
    v_rast <- tmp_spp_range_rast * v
    s_rast <- str_stack[[s]]
    i_rast <- v_rast * s_rast
    impact_list[[j]] <- i_rast
  }
  impact_stack <- stack(impact_list)
  impact_stack[['chi']] <- calc(impact_stack, sum)
  
  plot(impact_stack)

}

3.3.2 Dermochelys coriacea (leatherback turtle)

3.3.3 Carcharodon carcharias (great white shark)

3.3.4 Megaptera novaeangliae (humpback whale)

3.3.5 Thunnus orientalis (Pacific bluefin tuna)

3.3.6 Clupea harengus (Atlantic herring)

spp_vec <- spp_vuln_am  %>%
  filter(stressor %in% str_maps$str_vuln) %>%
  .$am_sid %>%
  unique()

spp_vuln_df <- spp_vuln_am %>%
  filter(am_sid %in% spp_vec) %>%
  select(am_sid, sciname, stressor, score, sd_score) %>%
  filter(stressor %in% str_maps$str_vuln)

spp_maps_df <- get_am_spp_cells(prob_cut = 0) %>%
  filter(am_sid %in% spp_vec) %>%
  mutate(presence = 1)
spp_chi_fstem <- here_anx('spp_impacts_hcaf/%s_chi_hcaf.csv')
### unlink(here_anx('spp_impacts_hcaf/*.*'))

n_done <- list.files(here_anx('spp_impacts_hcaf'), pattern = 'chi_hcaf.csv') %>%
  n_distinct()
tmp <- parallel::mclapply(spp_vec, mc.cores = 12,
    FUN = function(spp) {

      spp_chi_file <- sprintf(spp_chi_fstem, spp)

      if(file.exists(spp_chi_file)) {
        # message('File ', basename(spp_chi_file), ' exists... skipping!')
      } else {
        message('Processing impacts for ', spp)
        tmp_maps_df <- spp_maps_df %>%
          filter(am_sid == spp)
          
        tmp_spp_range_rast <- tmp_maps_df %>%
          map_to_hcaf(which = 'presence')
        tmp_spp_vuln_df <- spp_vuln_df %>%
          filter(am_sid == spp) %>%
          filter(score > 0)
        
        impact_list <- vector('list', length = nrow(tmp_spp_vuln_df)) %>%
          setNames(tmp_spp_vuln_df$stressor)
        for(j in 1:nrow(tmp_spp_vuln_df)) {
          ### j <- 1
          s <- tmp_spp_vuln_df$stressor[j]
          v <- tmp_spp_vuln_df$score[j]
          v_rast <- tmp_spp_range_rast * v
          s_rast <- str_stack[[s]]
          i_rast <- v_rast * s_rast
          impact_list[[j]] <- i_rast
        }
        impact_stack <- stack(impact_list)
        impact_stack[['chi']] <- calc(impact_stack, sum)
        
        
        impact_df <- as.data.frame(values(impact_stack)) %>%
          mutate(across(everything(), ~round(.x, 5))) %>%
          mutate(loiczid = 1:n()) %>%
          filter(!is.na(chi)) %>% 
            ### this works because all "present" cells have a numeric value (incl 0)
          left_join(tmp_maps_df, by = 'loiczid') %>%
          select(-am_sid, -presence)
        
        write_csv(impact_df, spp_chi_file)
      }
    })

3.4 Calculate ocean-area-weighted mean CHI for each species…

a_df <- get_hcaf_info() %>%
  select(loiczid, ocean_area) %>%
  distinct()

spp_chi_fs <- list.files(here_anx('spp_impacts_hcaf'), pattern = 'chi_hcaf.csv', full.names = TRUE)
spp_ids <- str_remove(basename(spp_chi_fs), '_chi_hcaf.csv')

mean_chi_list <- parallel::mclapply(spp_chi_fs, mc.cores = 16,
    FUN = function(f) {
      ### f <- spp_chi_fs[1]
      spp <- str_remove(basename(f), '_chi_hcaf.csv')
      spp_chi <- data.table::fread(f) %>%
        mutate(am_sid = spp) %>%
        left_join(a_df, by = 'loiczid') %>%
        group_by(am_sid) %>%
        summarize(area_km2 = sum(ocean_area),
                  mean_chi = sum(chi * ocean_area) / area_km2)
      if(nrow(spp_chi) == 0) {
        spp_chi <- data.frame(area_km2 = NA, mean_chi = NA)
      }
    }) %>%
  setNames(spp_ids)
mean_chi_df <- bind_rows(mean_chi_list, .id = 'am_sid')

write_csv(mean_chi_df, here('int/mean_chi_check.csv'))
mean_chi_df <- read_csv(here('int/mean_chi_check.csv')) %>%
  left_join(spp_vuln_by_gp %>% select(am_sid, p_group) %>% distinct(), by = 'am_sid')

ggplot(mean_chi_df, aes(x = log10(area_km2), y = mean_chi, color = p_group)) +
  geom_point(show.legend = FALSE) +
  geom_hline(yintercept = 0) +
  facet_wrap( ~ p_group) +
  theme_minimal()